DALS012-Linear Models线性模型02矩阵知识

前言

这一部分是《Data Analysis for the life sciences》的第5章线性模型的第2小节。这一部分的主要内容涉及矩阵的数学计算原理。

先来看一下前面的一个案例。

小鼠饲料案例

在前面的内容中，我们使用了t检验来研究高脂饲料(hf)饲喂的小鼠与普通饲料(chow)饲喂的小鼠体重的差异，现在我们使用线性模型来研究一下这两种小鼠的体重是否有差异，最终我们会发现，这两种统计方法在本质上是一样的，如下所示：

dir <- system.file(package = "dagdata")
list.files(dir)
dir
filename <- file.path(dir,"extdata/femaleMiceWeights.csv")
dat <- read.csv(filename)
# input raw data
head(dat)
str(dat)
stripchart(dat$Bodyweight ~ dat$Diet, vertical=TRUE, method="jitter",
main="Bodyweight over Diet")

结果如下所示：

> head(dat)
  Diet Bodyweight
1 chow      21.51
2 chow      28.14
3 chow      24.04
4 chow      23.45
5 chow      23.68
6 chow      19.79
> str(dat)
'data.frame':	24 obs. of  2 variables:
 $ Diet      : Factor w/ 2 levels "chow","hf": 1 1 1 1 1 1 1 1 1 1 ...
 $ Bodyweight: num  21.5 28.1 24 23.4 23.7 ...

从这张图上我们大致可以看出来，hf组的小鼠体重比chow组的小鼠体重高一些。现在我们使用设计矩阵$\mathbf{X}$（公式为~Diet）来计算一下这个结果，在设计矩阵中，第2列是分组信息，也就是饲料的类型，如下所示：

levels(dat$Diet)
X <- model.matrix(~ Diet, data=dat)
X

结果如下所示：

> X
   (Intercept) Diethf
1            1      0
2            1      0
3            1      0
4            1      0
5            1      0
6            1      0
7            1      0
8            1      0
9            1      0
10           1      0
11           1      0
12           1      0
13           1      1
14           1      1
15           1      1
16           1      1
17           1      1
18           1      1
19           1      1
20           1      1
21           1      1
22           1      1
23           1      1
24           1      1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$Diet
[1] "contr.treatment"

lm()函数背后的数学原理

现在我们先不用lm()来进行计算，先用设计矩阵%\mathbf{X}$来计算一下$\beta$，这个值能够使线性模型的平方和最小，公式如下所示：

$$\hat{\boldsymbol{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{Y}$$

在R中，可以使用矩阵相乘的符号%*%，以及solve()函数来求解上述方程，如下所示：

Y <- dat$Bodyweight
X <- model.matrix(~ Diet, data=dat)
solve(t(X) %*% X) %*% t(X) %*% Y

结果如下所示：

> Y <- dat$Bodyweight
> X <- model.matrix(~ Diet, data=dat)
> solve(t(X) %*% X) %*% t(X) %*% Y
                 [,1]
(Intercept) 23.813333
Diethf       3.020833

这些系数是对照组的平均值（23.813333），以及两组的差值（3.020833），计算一下，如下所示：

s <- split(dat$Bodyweight, dat$Diet)
mean(s[["chow"]])
mean(s[["hf"]]) - mean(s[["chow"]])

结果如下所示：

> s <- split(dat$Bodyweight, dat$Diet)
> mean(s[["chow"]])
[1] 23.81333
> mean(s[["hf"]]) - mean(s[["chow"]])
[1] 3.020833

现在我们使用lm()函数两周来计算一下，如下所示：

fit <- lm(Bodyweight ~ Diet, data=dat)
summary(fit)
(coefs <- coef(fit))

结果如下所示：

> summary(fit)

Call:
lm(formula = Bodyweight ~ Diet, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.1042 -2.4358 -0.4138  2.8335  7.1858 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   23.813      1.039  22.912   <2e-16 ***
Diethf         3.021      1.470   2.055   0.0519 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.6 on 22 degrees of freedom
Multiple R-squared:  0.1611,	Adjusted R-squared:  0.1229 
F-statistic: 4.224 on 1 and 22 DF,  p-value: 0.05192
> (coefs <- coef(fit))
(Intercept)      Diethf 
  23.813333    3.020833

线性回归系数

下图说明了前面我们计算出来的2个系数，即对照组小鼠的平均体重，以及hf组的小鼠体重与对照组的差异，如下所示：

stripchart(dat$Bodyweight ~ dat$Diet, vertical=TRUE, method="jitter",
           main="Bodyweight over Diet", ylim=c(0,40), xlim=c(0,3))
a <- -0.25
lgth <- .1
library(RColorBrewer)
cols <- brewer.pal(3,"Dark2")
abline(h=0)
arrows(1+a,0,1+a,coefs[1],lwd=3,col=cols[1],length=lgth)
abline(h=coefs[1],col=cols[1])
arrows(2+a,coefs[1],2+a,coefs[1]+coefs[2],lwd=3,col=cols[2],length=lgth)
abline(h=coefs[1]+coefs[2],col=cols[2])
legend("right",names(coefs),fill=cols,cex=.75,bg="white")

在这里使用前面的小鼠案例主要是为了说明了回归与t检验在本质上是一样的，这个简单的线性回归模型给出了与特定类型t检验相同的结果（t统计量与p值）。这是两组之间的t检验，前提是两组的总体标准差相同。当我们假设误差$\boldsymbol{\varepsilon}$都是均匀分布时，这些误差被编码到线性模型中。虽然在这种情况下的线性模相当于t检验，但我们可以对线性模型进行拓展，将其扩展到复杂的形式。在下面的内容中，我们将会说明线性模型能够得到几乎相同的结果：

> summary(fit)$coefficients
             Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 23.813333   1.039353 22.911684 7.642256e-17
Diethf       3.020833   1.469867  2.055174 5.192480e-02

t检验的结果与之相同：

ttest <- t.test(s[["hf"]], s[["chow"]], var.equal=TRUE)
summary(fit)$coefficients[2,3]
ttest$statistic

结果如下所示：

> summary(fit)$coefficients[2,3]
[1] 2.055174
> ttest$statistic
       t 
2.055174

练习

P189

标准误

使用矩阵代数来计算最小二乘估计时，所得到的估计值是取机变量。我们还能计算这些值的标准误。线性代数也能满足这些任务，先看一些案例。

自由落体

在学习统计学的过程中，了解随机源于何处是有必要的。在自由落体这个案例中，当我们对落下来的球进行测量时，就会出现检测误差，这就是自由落体运行中随机的来源。假设我们做了好几次实验，每次实验我们都会得到一组测量误差。这就意味着，我们得到的数据基本是随机变化的，这反过来，也意味着，我们计算得到的估计值也是随机变化的。在这个案例中，每次我们进行实验时，引力常数的估计值也都在发生变化。引力常数是一个确定的值，但是我们对它的估计不是。为了说明这一步，我们可以使用Monte Carlo模拟一下。具体来说，我们将会重复地生成数据，并对每次的二次项(quadratic term)进行估计，如下所示：

set.seed(1)
B <- 10000
h0 <- 56.67
v0 <- 0
g <- 9.8 ##meters per second
n <- 25
tt <- seq(0,3.4,len=n) ##time in secs, t is a base function
X <-cbind(1,tt,tt^2)
##create X'X^-1 X'
A <- solve(crossprod(X)) %*% t(X)
betahat<-replicate(B,{
y <- h0 + v0*tt - 0.5*g*tt^2 + rnorm(n,sd=1)
betahats <- A%*%y
return(betahats[3])
})
head(betahat)

结果如下所示：

> head(betahat)
[1] -5.038646 -4.894362 -5.143756 -5.220960 -5.063322 -4.777521

从上面结果我们可以发现，每次的估计值都不一样，这是因为估计值$\hat{\beta}$是一个随机变量，它其实是服从正态分布的，如下所示：

library(rafalib)
mypar(1,2)
hist(betahat)
qqnorm(betahat)
qqline(betahat)

上图显示的是，通过Monte Carlo模拟生成的自由落体运行数据对回归系数的估计值的分布。左侧是直方图，右侧是qq图。

由于$\hat{\beta}$是我们生成的那些服从正态分布的数据的线性组合，因此从QQ图上我们也可以看出，$\hat{\beta}$也服从正态分布。此外，这个分布的均值的真实参数0.5g，可以通过Monte Carlo模拟来所证实，如下所示：

round(mean(betahat),1)

结果如下所示：

> round(mean(betahat),1)
[1] -4.9

但是，当我们对估计值进行计算时，无法得到精确的数值，这是因为我们估计值的标准误大约是：

> sd(betahat)
[1] 0.2129976

在接下来的案例中，我们将会演示一下，不通过Monte Carlo模拟来计算标准误的方法，因为在实际的分析中，我们无法精确地知道误差是如何产生的。

父子身高

在父母身高的案例中，我们也遇到了随机问题，因为我们是对父子身高的总体进行随机抽样。为了说明这个问题，现在假设我们已经有了这个总体：

library(UsingR)
x <- father.son$fheight
y <- father.son$sheight
n <- length(y)

现在我们使用Monte Carlo模拟来生成一个含有50个数据的样本，如下所示：

N <- 50
B <-1000
betahat <- replicate(B,{
index <- sample(n,N)
sampledat <- father.son[index,]
x <- sampledat$fheight
y <- sampledat$sheight
lm(y~x)$coef
})
betahat <- t(betahat) #have estimates in two columns

现在绘制一个QQ图，我们看到，我们的估计值大致是服从标准正态分布的，如下所示：

mypar(1,2)
qqnorm(betahat[,1])
qqline(betahat[,1])
qqnorm(betahat[,2])
qqline(betahat[,2])

现在看一下估计值之间的相关性，如下所示：

> cor(betahat[,1],betahat[,2])
[1] -0.9992293

当我们计算我们估计值的线性组合时，我们需要知道这个信息能够正确地计算出这些线性组合的标准误差。在接下来的部分中，我们会提到方差-协方差(variance-covariance)矩阵。两个随机变量的协方差定义如下：

mean( (betahat[,1]-mean(betahat[,1] ))* (betahat[,2]-mean(betahat[,2])))

结果为：

> mean( (betahat[,1]-mean(betahat[,1] ))* (betahat[,2]-mean(betahat[,2])))
[1] -1.035291

协方差是相关系数乘以每个随机变量的标准差，如下所示：

$$\mbox{Corr}(X,Y) = \frac{\mbox{Cov}(X,Y)}{\sigma_X \sigma_Y}$$

除此之外，这个量在实际分析中没有一个现实意义上解释。但是，它对于数学推导的过程来说，有着重要意义。在接下来的部分中，我们会描述用于计算线性模型估计的标准差的矩阵代数计算方法。

方差-协方差矩阵

我们先来定义一什么是方差-协方差矩阵(variance-covariance matrix)$\Sigma$。地于一个元素是随机变量的向量$\mathbf{Y}$来说，我们定义$\Sigma$是含有$i,j$维的矩阵，如下所示：

$$\Sigma_{i,j} \equiv \mbox{Cov}(Y_i, Y_j)$$

如果$i=j$，那么协方差就等于方差，如果变量都是独立的，则协方差都为0。在我们现在所学到的各种向量里，$\mathbf{Y}$向量的每个观测值$Y_{i}$是从总体中抽样得到的，我们可以假设它的每个元素都是独立的，并且$Y_{i}$有着相同的方差$\sigma^2$，因此方差-协方差矩阵只有两类元素，如下所示：

$$\mbox{Cov}(Y_i, Y_i) = \mbox{var}(Y_i) = \sigma^2$$ $$\mbox{Cov}(Y_i, Y_j) = 0, \mbox{ for } i \neq j$$

这就是暗示了$\boldsymbol{\Sigma} = \sigma^2 \mathbf{I}$ ，其中$\mathbf{I}$是单位矩阵(Identity matrix)。

线性组合的方差

线性代数提供的一个有用结果就是$\mathbf{Y}$的$\mathbf{AY}$线性结合的方差-协方差矩阵，它的计算如下所示：

$$\mbox{var}(\mathbf{AY}) = \mathbf{A}\mbox{var}(\mathbf{Y}) \mathbf{A}^\top$$

例如，如果$Y_{1}$和$Y_{2}$是独立的，并且其方差为$\sigma^2$，那么：

$$\mbox{var}\{Y_1+Y_2\} = \mbox{var}\left\{ \begin{pmatrix}1&1\end{pmatrix}\begin{pmatrix} Y_1\\Y_2\\ \end{pmatrix}\right\}=\begin{pmatrix}1&1\end{pmatrix} \sigma^2 \mathbf{I}\begin{pmatrix} 1\\1\\ \end{pmatrix}=2\sigma^2$$

我们会使用上述的结果来获取LSE的标准误。

LES标准误

$\hat{\beta}$是一个线性组合，即 $\mathbf{Y}$: $\mathbf{AY}$ with $\mathbf{A}=\mathbf{(X^\top X)^{-1}X}^\top$的线性组合，因此我们可以使用上面的公式来计算一个我们估计值的方差：

$$\mbox{var}(\boldsymbol{\hat{\beta}}) = \mbox{var}( \mathbf{(X^\top X)^{-1}X^\top Y} ) = \\\mathbf{(X^\top X)^{-1} X^\top} \mbox{var}(Y) (\mathbf{(X^\top X)^{-1} X^\top})^\top = \\ \mathbf{(X^\top X)^{-1} X^\top} \sigma^2 \mathbf{I} (\mathbf{(X^\top X)^{-1} X^\top})^\top = \\ \sigma^2 \mathbf{(X^\top X)^{-1} X^\top}\mathbf{X} \mathbf{(X^\top X)^{-1}} = \\ \sigma^2\mathbf{(X^\top X)^{-1}}$$

其中，这个矩阵的平方根的对角线含有我们估计值的标准误。

估计$\sigma^2$

为了从上述公式中获得一个精确的估计值，我们需要估计$\sigma^2$。在前面我们通过抽样估计了标准误。但是$Y$的抽样标准误不是$\sigma$，因为$Y$还包含了变异，这些变量是由模型$\mathbf{X\beta}$的确定部分引入的。此时我们采用的方法是利用残差。

我们可以按下面的公式构建残差：

$$\mathbf{r}\equiv\boldsymbol{\hat{\varepsilon}} = \mathbf{Y}-\mathbf{X}\boldsymbol{\hat{\beta}}$$

其中，$r$和$\hat{\epsilon}$都可以用来表示残差(residuals)。然后我们使用上述的公式来估计残差，这种方式与单变量情况下类似：

$$s^2 \equiv \hat{\sigma}^2 = \frac{1}{N-p}\mathbf{r}^\top\mathbf{r} = \frac{1}{N-p}\sum_{i=1}^N r_i^2$$

其中，$N$是样本数目，$p$是$\mathbf{X}$的列数，或者是参数的数目（包括截矩$\beta_{0}$）。公式中还要除以$N-p$，这是因为根据数学理论（不用深究，这个跟自由度能关），除以$N-p$，表示无偏估计。下面在R中计算一下，观察一下我们是否能够获得与Monte Carlo模拟一样的数值：

n <- nrow(father.son)
N <- 50
index <- sample(n,N)
sampledat <- father.son[index,]
x <- sampledat$fheight
y <- sampledat$sheight
X <- model.matrix(~x)
N <- nrow(X)
p <- ncol(X)
XtXinv <- solve(crossprod(X))
resid <- y - X %*% XtXinv %*% crossprod(X,y)
s <- sqrt( sum(resid^2)/(N-p))
ses <- sqrt(diag(XtXinv))*s
summary(lm(y~x))$coef[,2]
ses
apply(betahat,2,sd)

结果如下所示：

> summary(lm(y~x))$coef[,2]
(Intercept)           x 
  8.3899781   0.1240767 
> ses
(Intercept)           x 
  8.3899781   0.1240767 
> apply(betahat,2,sd)
(Intercept)           x 
  8.3817556   0.1237362

从上面我们可以看出来，我们的计算结果与Monte Carlo模拟的结果几乎一样。

估计值的线性组合

多数情况下，我们会计算估计值的一个线性组合，例如$\hat{\beta_{2}}-\hat{\beta_{1}}$。这就是$\hat{\beta}$的一个线性组合，如下所示：

$$\hat{\beta}_2 - \hat{\beta}_1 = \begin{pmatrix}0&-1&1&0&\dots&0\end{pmatrix} \begin{pmatrix} \hat{\beta}_0\\ \hat{\beta}_1 \\ \hat{\beta}_2 \\ \vdots\\ \hat{\beta}_p \end{pmatrix}$$

通过上述的公式，我们就知道中如何计算$\hat{\beta}$的方差-协方差矩阵。

CLT与t分布

我们已经了解了如何过计算估计值的标准误。但是，我们在第1章中了解到，在计算置信区间时，我们需要知道随机变量的分布。我们努力计算标准误的原因是因为，CLT适用于线性模型。如果$N$足够大，那么LSE就会服从正态分布，其中这个正态分布的均值为$\beta$，标准误就是前面计算的标准误。对于小样本来说，如果$\varepsilon$服从正态分布，那么$\hat{\beta}-\hat{\beta}$服从t分布。在这里，我们不需要知道空上计算过程，但是这个结果很有，因为这是我们在线性模型中计算p值，置信区间的基础。

代码与数学

构建线性模型的标准做法要么是假设$\mathbf{X}$是固定的，要么我们给它们设置条件。因此$\mathbf{X\beta}$没有像$\mathbf{X}$那样固定的方差。这就浊为什么我们要写 $\mbox{var}(Y_i) = \mbox{var}(\varepsilon_i)=\sigma^2$。在实际计算中，这有可能造成误解，如下所示：

x = father.son$fheight
beta = c(34,0.5)
var(beta[1]+beta[2]*x)

结果如下所示：

> var(beta[1]+beta[2]*x)
[1] 1.883576

这个值并不等于0，也不接近于0。这就是为什么我们要谨慎区分代码与数学的一个案例。var()函数仅仅是用于计算我们输入数值的方差，而数学对方差的定义则是要考虑到这些数字是否是随机变量。在上述的R代码中，x并没有固定：我们只是让它发生变化，但是，当我们写下$\mbox{var}(Y_i) = \sigma^2$时，我们是把数学上的x强制固定下来。同样，如果我们使用R来计算自由落体运行案例中$Y$的方差，我们也会获得一个与$\sigma^2=1$（已知方差）明显不同的数值，如下所示：

n <- length(tt)
y <- h0 + v0*tt - 0.5*g*tt^2 + rnorm(n,sd=1)
var(y)

结果如下所示：

> n <- length(tt)
> y <- h0 + v0*tt - 0.5*g*tt^2 + rnorm(n,sd=1)
> var(y)
[1] 329.5136

练习

P199